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I. INTRODUCTION 



In recent papers Karmarkar [Refs. 1,2] has presented a new 
method for the solution of linear programming (LP) problems; His new 
solution technique moves from some feasible starting point across the 
interior region of a polytope that is defined by the problem constraints. 
He shows that the number of steps to find an optimal solution with his 
technique is polynomially bounded. 

In contrast, the simplex algorithm which, is widely used for the 
solution of LP problems finds an optimal solution by moving from vertex 
to vertex on the poly tope. This is known, in the worst case, to require 
an exponential number of steps [Ref. 3] . 

The polynomial bound of the projective algorithm makes the new 
solution technique very appealing to researchers. The theory of the 
new algorithm seems to be widely accepted among experts, while 
Karmarkar's claim that his algorithm is 50-100 times faster than the 
simplex method has met skepticism [Ref. 4]. 

In this study a variant of the projective method is implemented, 
and some well known test problems are solved. 

A. THE BASIC PROJECTIVE METHOD 

Following Karmarkar [Ref. 1: p. 4], the number of steps of the 
algorithm depends on R/r, where R is the radius of the sphere 
circumscribing the polytope, and r the radius of the inscribed sphere. 

Assume a general LP of the form 

T 

Min c x x 

s.t. A x = b (LP1) 

x > 0 . 

With the assumption that the sum of its variables has an upper bound, 
and with the proper scaling of variables, a convexity constraint 
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( 1 . 1 ) 




can be added to LP1. This transformation maps the LP onto a unit 
simplex S whose center is at a Q = (1/n, . . . , 1/n) . Then, a transforma- 
tion is performed such that a feasible interior starting point is mapped 
onto a Q . 

If B(a Q ,r) is the largest sphere with center a Q that can be 
inscribed into the simplex S, and -B(a o ,R) is the smallest sphere 
circumscribing S, then 



By restricting a solution x to remain inside the largest inscribed sphere 
B(a Q ,r), the method achieves in one iteration a reduction in the differ- 
ence between the current objective value and the optimal objective value 
by a factor of (1 - 1/n). Following Shanno [Ref. 5] a simple proof is: 



R/r = n. 



( 1 . 2 ) 



Let 



f' - min c T x, x 6 S, Ax = b, 



(1.3) 



T T 

X = min c x, x e B(a Q ,r), Ax = b, l^x = 1, 



(1.4) 



f = min c T x, x £ B(a Q ,R), Ax = b, l T x = 1. 



(1.5) 



Then 



c T a Q - _f < c T a 0 - f" < c T a Q - f, 



( 1 . 6 ) 



and with equation 1.2 



c T a 0 - f = n(c T a 0 - X) . 



(1.7) 



From that 



(X - O/(c T a 0 - O < (1 - 1/n). 



( 1 . 8 ) 
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If a linear objective function could be maintained at each iteration, 
it follows that an upper bound on the number of steps required to find 
an optimal solution is of 0(n). Unfortunately, the projective transfor- 
mations needed to continue the algorithm result in a nonlinear objective 
function. 

B. DERIVATION OF THE ALGORITHM 
1. The Canonical Form 

Suppose we have a linear programming problem of the form 
LP1. Karmarkar*s method requires that this LP be transformed into the 
following canonical form, 



= 0 (LP2) 

= 1 
> 0 , 

where the optimal solution value is zero. With the assumption that the 
sum of the variables of an LP is bounded above and subsequent scaling 
by this bound, a convexity constraint (equation 1.1) can be added to 
LP1. In practice this can be a problem because choosing too large an 
upper bound may cause numerical problems. 

LP2 requires that the nonhomogeneous system of equations 
Ax = b be made homogeneous. Karmarkar [Ref. 1: p.34] proposes a 

transformation that would transform the i-th equation a* x = b^ to 

pa.j - b.)x. = 0. (1.9) 

This transformation has the disadvantage that when b is dense the 
sparsity of A will be lost. Karmarkar requires the optimal solution 
value to LP2 to be zero together with a special stopping criterion 
(equation 1.26), to prove the polynomial bound. To achieve the zero 
objective value the optimal objective function value f‘ of the 



Min 



T 

c x x 



s.t. Ax 
1 T X 
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untransformed LP has to be known in advance, and the following trans- 
formation made 



c T x - f* = C T X 



ri T x = (c -f*i) T x . 



( 1 . 10 ) 



Karmarkar concludes [Ref. 2: p. 387] that if the minimal objective value 
determined by the projective algorithm is not equal to zero, the original 
problem must be either infeasible or unbounded. Another method for 
making the transformation from LP1 to LP2 is discussed in 
Chapter II. B. 

2 . The Projective Transformation 

Let x° > 0 be a known feasible solution to LP2. The invertible 



projective transformation 



y = 



D _1 x 



1 T D _1 x 



( 1 . 11 ) 



l T x = 



1 and x > 0, onto y such .that 

is a diagonal matrix whose entries are 



maps any x, such that 
l^y = 1 and y > 0. D 
( Xl 0 ,...,x n 0 ). 

The transformation maps the LP2 unit simplex in x-space onto 



another unit simplex in y-space. The point 



is mapped onto 



y° = 1/n l 1 , the center of the unit simplex in y-space. The inverse of 
the transformation (equation 1.11) is given by 



X = 



D Y 
1 T D y 



( 1 . 12 ) 



After the transformation, LP2 can now be restated as 
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c T D y 

Min 

1 T D y 

s.t. A D y =0 (LP3) 

1 T y =1 
y > 0 . 

Strictly speaking LP3 is not a linear program because its objective 
function is a rational function of y. Karmarkar [Ref. 1: page 18ff] 

shows that it is sufficient to consider a linear approximation to the 
objective function of LP3, which leads to 

Min c^D y 

s.t. A D y =0 (LP4) 

1 T y =1 
y > 0 . 



3. Optimization Over a Sphere 

A solution to LP4 is now restricted to lie within a sphere with 
center at y° and radius «r, where 

r = l/(n(n-l))' 1/2 . (1.13) 

r is the radius of the largest sphere that can be inscribed into the unit 
simplex, and a is a constant such that 0<a<l. a provides a margin 
which ensures that the algorithm doesn't select a point outside the 
sphere due to round-off error. By convexity, an optimal solution will 
occur at the boundary of the sphere. Thus, the additional constraint 
can be stated as an equality. Now we have ^ the following version of the 
LP: 
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0 

1 



(LP5) 



c T D y 



Min 
s.t. A D y 

i T y 

(y-y°) T (y-y°) 



.2 2 

<X r 



Before the problem can be solved, one more transformation 
that moves the center of the sphere to the origin is useful. Let 
Y = y + y°, and eliminate the constant terms ADy° = 0, l^y 0 = 1 and 

m n 

c x Dy . For convenience, define B by 



A D 




and the gradient c of the objective function of LP5 by 
c = c^D . 



(1.14) 



(1.15) 



Then LP5 can be restated as, 



— T— 

Min c y 

s.t. By = 0 

— T— 2 2 

y 1 y = a z r z . 



(LP6) 



In order to solve LP6, note that an initial feasible solution to 
LP6 is y = 0 (which corresponds to y = y° in LP5). This is also the 

T 9 2 

center of the sphere y y = a r . An optimal solution to LP6 can be 
obtained by finding a direction of maximum rate of ascent c that is 
feasible with respect to By = 0, and moving in direction -c (maximum 
rate of descent) a distance from y = 0 to the boundary of the 

sphere . 

A feasible direction of maximum ascent is found by orthogo- 
nally projecting c onto the null space of B (see [Ref. 1: page 17] ), 
i.e., the following problem has to be solved in terms of Cp 
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(1.16) 



Min c T c - 2c T c + c p Tc p = Min ( || c - c p || 2 ) 2 

C P C P 

s.t. Bc p = ® • 

Define the Lagrangian L(c p , \) by 

L(c p ,70 = c T c - 2c T c p ♦ c p T c p - X T Bc p . (1.17) 

First-order optimality conditions for the Lagrangian yield 

-2c + 2c p - B T A = 0, (1.18) 

and 

Bc p = 0. (1.19) 

After multiplying by B and dropping 2Bc p =0 (equation 1.19) we get 

-2Bc = BB T X . (1.20) 

T - 1 

Assuming (BB ) exists, we can solve for X 



A = -2(BB T ) -1 Bc. (1.21) 

Substituting equation 1.21 into equation 1.18 gives 

C p = c - B T (BB T )' 1 Bc. (1.22) 

The direction of maximum rate of ascent is then 

C = C p/H CpH • (1.23) 

Thus, the optimal solutions to LP6 and LP5 are 

y = - «rc, (1.24) 
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and 



y 1 = y° - d re (1-25) 

respectively. 

Finding Cp is the key part of Karmarkar's method because it 
involves the major portion of the computational work. Solving equation 
1.22 for Cp can be viewed as solving a linear least- squares problem 
(see Chapter II. C.). 

With the optimal solution to LP5 in y- space, the method 
proceeds by transforming that solution back into x- space by use of 
equation 1.12. Next it defines a new matrix D and iterates until it 
reaches the stopping criterion [Ref. 1: p. 14] 

(c T x k )/(c T x°) < 2"^ (1.26) 

T 

where q is a termination parameter. Note that c x = 0 at optimality. 
Table 1 shows an outline of Karmarkar's basic method. 
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TABLE 1 










ALGORITHM 


1 




Input 


Problem Size 








. . .n 




Coefficient 


Matrix . . 




. . .A 




Cost 


Function . 






...c T 




Initial Feasible Solution 


...x° 




Termination 


Parameter. . . . 


. . .q 




Feasibility 


Parameter. . . . 


. . . c* 


Begin 


k = 


1 












x = 


x° 












f° = 


C T X° 












f = 


OO 












r = 


l/(n(n 


-1))' 1/2 




While 


(f/f° 


> 2‘3) 














D = 


diag (x) 








y = 


(D' 1 x)/(1 T D' 


1 x) = 1 






c = 


c T D 










B = 




’AD' 







C P 


= 


c - B T ( BB T ) " 1 Bc 


A 

C 


= • 


c p/ll c pll 


y 


r= 


y - arc 


X 


- 


(Dy )/ ( l T Dy ) 


f 


= 


T 

c A x 


k 


- 


k + 1 


End (While) 






End 






Output Solution.. 




X 


Obj ective 


Function Value... f 


Iteration 


Count r k 
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II. A VARIANT OF THE PROJECTIVE METHOD 



A. INTRODUCTION 

The following sections outline a variant of the projective method 
that was proposed by Wood [Ref. 6]. It has a practical method of 
bringing an LP into the canonical form, and uses the non-linear objec- 
tive function of LP3 to find a gradient c f which corresponds to c of 
equation 1.16. Similar to the simplex method, the proposed algorithm 
uses a ratio test to determine a feasible step length. 

Other practical features of the proposed method include the relax- 
ation of the requirement to know the optimal objective value in advance, 
and exploitation of sparsity of A. The implementation is especially 
concerned with controlling fill-in during the solution of the normal 
equations . 

B. DERIVATION OF THE MODIFIED ALGORITHM 

Given an LP problem of the form LP1, we use a single artificial 
variable + ^ to attain initial feasibility: 

Min (c T ,M)(x,x n + 1 ) 

s.t. Ax-(Ax°-b)x n + 1 = b (LP7) 

x,x n+ i ^ 0 . 

where M is the cost for the artificial variable and x° is the initial solu- 
tion with 

(x°,x°n,i) = 1 T . (2.1) 

The transformation has added one more column to the problem. 

To get a homogeneous right-hand side in LP7 we introduce an 
additional variable, and apply the following projective transformation, 
whose inverse is given by 
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(x, x n + i) = (n+2)(x',x' n+1 )/x' n+2 , 



( 2 . 2 ) 



(x',x' n + i) = (n+2) (x,x n + 1 )/(l T (x,x n+1 )+l) (2.3) 

x'n +2 = (n + 2)/(l T (x,x n+1 ) + l) . (2.4) 

Let the artificial column be denoted by a, i.e. a = -(Ax°-b). The 
transformed problem can now be restated, with the exception of the 
objective function, in Karmarkar's canonical form, 

(c T ,M, 0) (x',x' n + 1 ,x' n + 2 ) 
x 'n + 2 

(A, a, -b) (x',x , n+1 ,x' n+2 ) = 0 (LP8) 

l T (x , ,x' n+1 ,x' n+2 ) = n+2 

(x',x' n+1 ,x' n+2 ) > 0. 



Min 
s . t . 



The above transformation adds yet another column to the problem, but 
has the advantage that it doesn't change the sparsity of A, as opposed 
to Karmarkar's proposal. Rather than projecting the LP problem onto a 
unit simplex, it is projected onto an (n + 2) -simplex. This is done to 
improve numerical stability. 

The projective transformation, equation 1.12, is applied to LP8 
which gives 



Min (c T ,M,0)D(y,y n + 1 ,y n+2 )/d n+2 y n+2 

s.t- (A,a,-b)D(y,y n+1 ,y n+2 ) = 0 (LP9) 

i T (y,y n+ i,y n+2 ) = n+2 

(y,y n+ 1 ,y n+2 ) ^ 0 • 

The gradient (compare with equation 1.16) of the objective function of 

o T 

LP9, evaluated at the initial feasible starting point y = 1 is propor- 
tional to 
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c' = (c 1 d 1 ,...,c n d n ,Md n+1 ,-c T x). 



(2.5) 



The step of normalizing Cp in Algorithm 1 is replaced by a ratio 
test. Let 

jmax = argmax{(Cp)j). (2.6) 

This allows the algorithm to make a step outside the inscribed sphere, 
but maintains feasibility by restricting a solution to lie inside the 
simplex. Then, the update of y becomes 

Y - 1 - ^ c p/ ( c p)jmax ' ( 2 * 7 ) 

where p is a parameter to maintain feasibility such that 0<p<l. Table 2 
shows the modified algorithm. 

C. THE LINEAR LEAST- SQUARES PROBLEM 

Computation of the projected gradient Cp during every iteration 
accounts for most of the computational workload in any algorithm based 
on Karmarkar's method. Solving equation 1.22 can be viewed as solving 
the following linear least-squares problem, 

Min ( lie' - B T AII 2 ) 2 (2.8) 

T 

where B is an (n + 2)xm matrix with m<(n+2) assumed. 

T 

If rank(B )=m, then the solution to (2.8) is given by the solution 
to the system of normal equations 

BB T A = Be' . (2.9) 

The projected gradient c is then the residual vector of the least- 

r 

squares problem (2.8), i.e., 

c p = c' - B T A . (2.10) 
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Input 



TABLE 2 
ALGORITHM 2 

Problem Size n 

Coefficient Matrix A 

T 

Cost Function c 

Initial Feasible Solution ... x° 

Termination Parameter t 

Feasibility Parameter P 



Begin 



While 



k 

x 



= 1 



o _ 



f 
(f° 



End (While ) 



x 

oo 

f(x) 

D 

f° 

y 

c' 

B 



P 

^max 

y 

x 

k 



t ) 



diag (x) 
f (x) 

(n+2) (D" 

V (y) 



1 x)/(1 T D" 1 x) 



A D 
1 T 

= c' - B T (BB T ) _1 Bc' 
= argmax(Cp)j 

= 1 - ?c p/ (C P } imax 

= ( n+2 ) ( Dy ) /( l^Dy ) 

= k + 1 



End 



Output 



Solution x 

Objective Function Value... f 
Iteration Count k 
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